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POLYNOMIALS IN ASYMPTOTICALLY FREE RANDOM 

MATRICES 

ROLAND SPEICHER 


Abstract. Recent work of Belinschi, Mai and Speicher resulted in a general 
algorithm to calculate the distribution of any selfadjoint polynomial in free 
variables. Since many classes of independent random matrices become asymp¬ 
totically free if the size of the matrices goes to infinity, this algorithm allows 
then also the calculation of the asymptotic eigenvalue distribution of polyno¬ 
mials in such independent random matrices. We will recall the main ideas of 
this approach and then also present its extension to the case of the Brown 
measure of non-selfadjoint polynomials. 


1. Introduction 

Recent work of Belinschi, Mai and Speicher [2] resulted in a general algorithm 
to calculate the distribution of any selfadjoint polynomial in free variables. Since 
many classes of independent random matrices become asymptotically free if the size 
of the matrices goes to infinity, this algorithm applies then also to the calculation of 
the asymptotic eigenvalue distribution of polynomials in such independent random 
matrices. Here we will, after first recalling the main ideas of this approach, also 
address the non-selfadjoint situation and present the results of work in progress of 
Belinschi, Sniady and Speicher [1]. There it is shown that the combination of the 
ideas from the selfadjoint case with the hermitization method allows to extend our 
algorithm also to the calculation of the Brown measure of arbitrary polynomials 
in free variables. It is expected that this Brown measure is then also the limiting 
eigenvalue distribution of corresponding random matrix problems; however, since 
the convergence of ^-moments does not necessarily imply the convergence of Brown 
measures in the non-normal case, this has to remain as a conjecture for the moment. 


2. The Case of One Matrix 

We are interested in the limiting behaviour of N x N random matrices for N —» 
oo. Let us first recall the well-known situation of one matrix. 

Typical phenomena for basic random matrix ensembles are that we have al¬ 
most sure convergence of the eigenvalue distribution to a deterministic limit and, 
furthermore, this limit distribution can be effectively calculated. 


Talk given at the conference “Random Matrix Theory: Foundations and Applications” in 
Krakow, Poland 2014. 

This work has been supported by the ERC Advanced Grant “Noncommutative Distributions in 
Free Probability” NCDFP 339760. 


1 



2 


R. SPEICHER 


The common analytic tool for calculating this limit distribution is the Cauchy 
transform. For any probability measure fionR we define its Cauchy transform 



R 


This is an analytic function G : C + —> C and we can recover /i from G by Stieltjes 
inversion formula 

d/j(t) = — — lim 3G(f + ie)dt. 

7T e->0 

Quite often, one prefers to work with —G, which is called the Stieltjes transform. 

For the basic matrix ensembles, the Gaussian and the Wishart random matri¬ 
ces, one can derive quadratic equations for the respective Cauchy transforms and 
thus get an explicit formula for G(z). In Figure[l]we show the comparision between 
histograms of eigenvalues and the theoretical limit result given by the Cauchy trans¬ 
form. 




Figure 1. left: histogram of the eigenvalues of a 4000 x 4000 
Wigner matrix, compared to Wigner’s semicircle, determined by 
its Cauchy transform G{z) = z ^ v '^ 2 ~ 4 ; right: histogram of the 
eigenvalues of a Wishart matrix XX* for a 3000 x 12000 matrix 
X, compared to the Marchenko-Pastur distribution for A = 1/4, 
determined by its Cauchy transform G(z) = 

2+1 — A— yj (2 —(1+A)) 2 —4A 
22 


3. The Case of Several Independent Matrices 

Instead of one random matrix we are now interested in the case of several in¬ 
dependent random matrices. Since there is no meaningful notion of joint eigen¬ 
value distribution of several non-commuting matrices, we are looking instead on 
the eigenvalue distribution of polynomials in those matrices. For the moment, we 
restrict to selfadjoint polynomials in selfadjoint matrices. Later we will look on the 
non-selfadjoint case. So we are interested in the limiting eigenvalue distribution 
of general selfadjoint polynomials p( A' 1; ..., Xyj) of several independent selfadjoint 
N x N random matrices Xi ,..., X 
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Again one observes the typical phenomena of almost sure convergence to a de¬ 
terministic limit eigenvalue distribution. However, this limit distribution can be 
effectively calculated only in simple situations. 

In Figure [2] we show two such typical situations. The first is the polynomial 
p(X, Y) = X + Y, for independent Wigner matrix X and Wishart matrix Y; in this 
case, one can derive quite easily the implicit equation G(z) = Gwishart(-z — G(z)) for 
the Cauchy transform G{z) of the limiting distribution oip(X, Y); this subordinated 
form of equation is quite amenable to numerical solutions via iterations and serves 
as the model what we can hope for in more general situations. The second case 
is the polynomial p(X,Y) = XY + YX + X 2 , for which no theoretical result has 
existed up to now. 



Figure 2. Histogram for a generic realization of a 3000 x 3000 
random matrix p(X, Y), where X and Y are independent Gaussian 
and, respectively, Wishart random matrices: p(X, Y) = X + Y 
(left); p(X,Y) = XY + YX + X 2 (right). In the left case, the 
asymptotic eigenvalue distribution is relatively easy to calculate; 
in the right case, no such solution was known, this case will be 
reconsidered in Figure [3] 

There existed a huge literature on special cases of such polynomials in indepen¬ 
dent matrices (see, for example, [51 E H H ES [13 EH), but up to |2] there has 
been no algorithm for addressing the general problem. 

4. Asymptotic Freeness of Random Matrices 

The problem of the asymptotic eigenvalue distribution of random matices became 
linked to free probability theory by the basic result of Voiculescu [20j that large 
classes of independent random matrices (like Wigner or Wishart matrices) become 
asymptoticially freely independent. 

The conclusion of this is: calculating the asymptotic eigenvalue distribution of 
polynomials in such matrices is the same as calculating the distribution of polyno¬ 
mials in free variables. The latter is an intrinsic problem within free probability 
theory. 
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We want to understand the distribution of polynomials in free variables. What 
we understand quite well, by the analytic theory of free convolution, is the sum of 
free selfadjoint variables. So we should reduce the problem of arbitrary polynomials 
to the problem of sums of selfadjoint variables. 

This sounds like a quite ambitious goal, but it can indeed be achieved. However, 
there is a price to pay: we have to go over to an operator-valued frame. 


5. The Operator-Valued Setting 

Let A be a unital algebra and B C A be a subalgebra containing the unit. A 
linear map E : A B is a conditional expectation if 

E[b] =b \/b£B 

and 

E[biab 2 \ = biE[a\b 2 Va £ A, Vbx,b 2 £ B. 

An operator-valued probability space consists of B C A and a conditional expec¬ 
tation E : A -A B. Then, random variables Xi £ A (i £ I) are free with respect 
to E (or free with amalgamation over B) if E[a\ ■ ■ ■ a„] = 0 whenever cq £ B(xj m) 
are polynomials in some 2y(q with coefficients from B and E[af\ = 0 for all i and 
j(!) j( 2 ) 7^ / j ( n ). 

For a random variable x £ A, we define the operator-valued Cauchy transform: 

G(b) := E[{b - x)- 1 ] (4eB), 

whenever ( b — x ) is invertible in B. 

In order to have some nice analytic behaviour, we will in the following assume 
that both A and B are C*-algebras; B will usually be of the form B = M N { C), the 
N x TV-matrices. In such a setting and for x = x*, this G is well-defined and a nice 
analytic map on the operator-valued upper halfplane: 

H +(B) := {b £ B | (6 - b*)/(2i) > 0} 


and it allows to give a nice description for the sum of two free selfadjoint elements. 
In the following we will use the notation 


h(b ) := 


1 

W) 


- b. 


Theorem 1 ([2]). Let x and y be selfadjoint operator-valued random variables free 
over B. Then there exists a Frechet. analytic map u>: H + (H) —> H + (H) so that 


G x+y (b) = G x (uj(b)) for all b £ H+(B). 


Moreover, if b £ then ui(b) is the unique fixed point of the map 


f b : m+(B)^m+(B), f b (w) = hy(h x (w) + b) + b , 


u>(b) = lim f b n (w ) for any w £ H + (S). 

n—>-oo 


and 
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6. The Linearization Philosophy 

As we just have seen, we can deal with the sum of free selfadjoint elements, even 
on the operator-valued level. However, what we are interested in are polynomials 
in free variables, on a scalar-valued level. The relation between these two problems 
is given by the linearization philosophy: in order to understand polynomials p in 
non-commuting variables, it suffices to understand matrices p of linear polynomials 
in those variables. 

In the context of free probability this idea can be traced back to the early papers 
of Voiculescu; it became very prominent and concrete in the seminal work [5] of 
Haagerup and Thorbjprnsen on the largest eigenvalue of polynomials in indepen¬ 
dent Gaussian random matrices. A more streamlined version of this, based on the 
notion of Schur complement, which also has the important additional feature that 
it preserves selfadjointness, is due to Anderson [ 1 ] 

Actually, those ideas were also used before in the context of automata theory 
and formal languages or non-commutative rational functions, where they go under 
different names, like descriptor realization. Some relevant literature in this context 
is ns la ns]. Actually, from this context it becomes clear that we cannot only 
linearize polynomials, but also non-commutative rational functions. 

The crucial point is that for any selfadjoint polynomial p there exists a lineariza¬ 
tion p, which is also selfadjoint. Since p is linear it is the sum of operator-valued 
variables. 

We will present the details of this procedure with the help of a concrete example. 
Let us consider the polynomial p(x, y ) = xy + yx + x 1 in the free variables x and y. 
This p has a linearization 

/ 0 x y+f\ 

p=\ x 0 1 > 

\y +1 -i o ) 

which means that the Cauchy transform of p can be recovered from the operator¬ 
valued Cauchy transform of p, namely we have 

G p (b) = id 0 <p(( 6 -p)“ 1 ) = ^ for b = o cfj . 

But this p can now be written as 

(0 x f\ (0 0 y\ 

p = lx 0 -1 + 0 0 0 \ = X + Y 

If -1 0/ \y o 0 J 

and hence is the sum of two selfadjoint variables X and Y , which are, by basic 
properties of freeness, free over C). So we can use our subordination result in 
order to calculate the Cauchy transform of p. 

(' ~ p)_1) ;;;) = g*(&) = c^b) = g *m&)), 

where we calculate u>(b) via iterations as in Theorem [l] Figure [3] shows the agree¬ 
ment between the achieved theoretic result and the histogram of eigenvalues from 

Fig. [3 
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Figure 3. Plots of the distribution oip{x,y) = xy + yx + x 2 (left) 
for free x, y , where x is semicircular and y Marchenko-Pastur, and 
of the rational function r(x\,x 2 ) (right) for free semicircular ele¬ 
ments xi and x 2 ; in both cases our theoretical limit curve is com¬ 
pared with the histogram of the corresponding eigenvalue problem; 
for the left case compare Figure [2] 


Conclusion: the combination of linearization and operator-valued subordination 
allows to deal with case of selfadjoint polynomials in free variables, and thus with 
selfadjoint polynomials in asymptotically free random matrices. 


7. The Case of Rational Functions 


As we mentioned before this linearization procedure works as well in the case 
of non-commutative rational functions (and this was actually the main object of 
interest in the work of Schiitzenberger). This is work in progress with Mai, let us 
here just give an example for this. 

Consider the following selfadjoint rational function 

r(x!,x 2 ) = (4 - xi) _1 + (4-x 1 )~ 1 x 2 ((4 - xx) - x 2 (4 - x 1 )~ 1 x 2 ) 1 x 2 (4 - xi) _1 


in two free variables X\ and x 2 . The fact that we can write it as 

\xi y 1 

-\x 2 \-\xx) 

gives us immediately a selfadjoint linearization of the form 

0 


r(x\,x 2 ):={\ 0) 


r(x i,x 2 ) = 


i 

4^2 


\x 2 


-1 + 4 a 'l > 


(0 

1 

9 

0 ^ 

, (o 

0 

0 \ 

u 

' + 

0 | 

+ [\ 

0 


\o 

0 

-1+ \xx) 

V° 

\X2 

0 ) 


In Figure[3j we compare the histogram of eigenvalues of r(X-|, X 2 ) for one realiza¬ 
tion of independent Gaussian random matrices Xx,X 2 of size 1000 x 1000 with the 
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distribution of r{x 1 , 22 ) for free semicircular elements X\,X 2 , calculated according 
to our algorithm. 


8. The Brown Measure for Non-Normal Operators 

The main point of this note is that we can actually also extend the previous 
method to deal with not necessarily self-adjoint polynomials in free variables. 

The first problem in this situation is that we need to describe a measure on the 
complex plane by some analytic object. Let us first have a look on this. 

For a measure /i on C its Cauchy transform 

G>( A) = [ y— dp(z) 

Jc z 

is well-defined everywhere outside a set of R 2 -Lebesgue measure zero, however, it 
is analytic only outside the support of /i. 

The measure f can be extracted from its Cauchy transform by the formula 
(understood in distributional sense) 

A better approach to this is by regularization: 

GeA X ) = [ 2 , A | x " 12 d v{ z ) 

J c e + I A --I 

is well-defined for every A € C. By sub-harmonicity arguments 


nd% G ^ {X) 

is a probability measure on the complex plane and one has the weak convergence 
lim £ ^o AL = M- 

Our general polynomial in free variables will in general not be selfadjoint nor 
normal. Thus we need a generalization of the distribution of a selfadjoint operator 
to a non-normal situation. This can be done, in a tracial setting, by imitating 
essentially the ideas from above. 

So let us consider a general (not necessarily normal) operator i in a tracial non- 
commutative probability space ( A , ip). (We need some nice analytic frame here, so 
A should be a von Neumann algebra.) We put 

G e , x ( A) := ip ((A - s)*((A - *)(A - x)* + e 2 )” 1 ) . 


Then 


f^e,x 


ld_ 

7T d\ 


GcM) 


is a probability measure on the complex plane, which converges weakly for e —> 0. 
This limit f x := lim^o is called the Brown measure of x ; it was introduced 
by L. Brown in 1986 and revived in 2000 by Haagerup and Larsen [7], who made 
decisive use of it for their investigations around the invariant subspace problem. 
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9. Hermitization Method 


The idea of the hermitization method is to treat non-normal operators (or ran¬ 
dom matrices) x by studying sufficiently many selfadjoint 2x2 matrices built out 
of x. A contact of this idea with the world of free probabilty was made on a formal 
level in the works of Janik, Nowak, Papp, Zahed (TT! and of Feinberg, Zee [6]. We 
show in [4] that operator-valued free probability is the right frame to deal with 
this rigorously. By combining this with our subordination formulation of operator¬ 
valued free convolution we can then calculate the Brown measure of any polynomial 
in free variables. 

In order to get the Brown measure of x, we need 

G e ,x(A) = ((A - x)*((A - a;)(A - a;)* + e 2 ) ^ . 


Let 

x = (" 'o) eA «A. 

Note that A' is selfadjoint. Consider A in the M 2 ((C)-valued probability space with 
repect to E = id (g) tp : M 2 (A) —> M 2 (C) given by 


{ CL 11 

ai 2 \ 

A a 2i 

a 2 2/ 


(<p(an) V?(ai 2 )\ 
V^(«2l) ^( 022 )/ 


For the argument 

A «=(x fe) eM ^ c ) 

consider now the M 2 (C)-valued Cauchy transform of X 


Gx(Ae)=E[(A e -X)~ 1 ] 


( ffe, A,11 ge, A,12 
\ffe,A,21 Se, A,22 


One can easily check that (A e — X) 1 is actually given by 


( — ie((X — x)(A — x)* + e 2 ) 1 (A — s)((A — x)*(X — x) + e 2 ) 

\^(A — x)*((A — a:)(A — a;)* + e 2 ) -1 —*e((A — a;)*(A — a;) + e 2 ) -1 J 


and thus we are again in the situation that our quantity of interest is actually one 
entry of an operator-valued Cauchy transform: 


ge, a, 12 = G £iX ( A). 


10. Calculation of the Brown Measure 

So in order to calculate the Brown measure of some polynomial p we should 
first hermitize the problem by going over to selfadjoint 2x2 matrices over our 
underlying space, then we should linearize the problem on this level and use finally 
our subordination description of operator-valued free convolution to deal with this 
linear problem. It might be not so clear whether hermitization and linearisation 
go together well, but this is indeed the case. Essentially we do here a linearization 
of an operator-valued model instead of a scalar-valued one: we have to linearize a 
polynomial in matrices. But the linearization algorithm works in this case as well. 

Let us illustrate this with an example. 

Consider the polynomial p = xy in the free selfadjoint variables x = x* and 
V = V*- 
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For the Brown measure of this we have to calculate the operator-valued Cauchy 
transform of 

0 xy^ 

K yx 0 

In order to linearize this we should first write it as a polynomial in matrices of 
x and matrices of y. This can be achieved as follows: 


P = 


P = 


0 xy 
yx 0 


x 0 
0 1 


0 y 

y 0 


x 0 
0 1 


= XYX. 


P = XYX is now a selfadjoint polynomial in the selfadjoint variables 


X =(S 5) “ d y =(“ ») 

and has thus a selfadjoint linearization 

/ 0 0 X\ 

b r 1 ■ 

U 1 o/ 

Pluggin in back the 2x2 matrices for X and Y we get finally the selfadjoint 
linearization of P as 


(0 

0 

0 

0 

X 

0 \ 

0 

0 

0 

0 

0 

1 

0 

0 

0 

y 

-1 

0 

0 

0 

y 

0 

0 

-1 

X 

0 

-l 

0 

0 

0 

V° 

1 

0 

-l 

0 

0 / 


which can be written as the sum of two Mq(C )-free matrices: 


(0 

0 

0 

0 

X 

0 ^ 


/o 

0 

0 

0 

0 

°\ 

0 

0 

0 

0 

0 

1 


0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

-1 

0 


0 

0 

0 

y 

0 

0 

0 

0 

0 

0 

0 

-1 

+ 

0 

0 

y 

0 

0 

0 

X 

0 

-1 

0 

0 

0 


0 

0 

0 

0 

0 

0 

v° 

1 

0 

-1 

0 

0 / 


v> 

0 

0 

0 

0 

0/ 


For calculating the Cauchy transform of this sum we can then use again our subordi¬ 
nation algorithm for the operator-valued free convolution from Theorem [l] Putting 
all the steps together gives an algorithm for calculating the Brown measure of p. 

Of course, we expect that in nice cases the eigenvalue distribution of our polyno¬ 
mial evaluated in independent Wigner or Wishart matrices should converge to the 
Brown measure of the polynomial in the corresponding free variables. However, in 
contrast to the selfadjoint case this is not automatic from the convergence of all 
relevant moments and one has to control probabilities of small eigenvalues during 
all the calculations. Such control have been achieved in special cases (in particular, 
the circular law and the single ring theorem), but in general it has to remain open 
for the moment. 

In the following figures we give for some polynomials the Brown measure accord¬ 
ing to our algorithm and compare this with histograms of the complex eigenvalues 
of the corresponding polynomials in independent random matrices. 
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As before, this algorithm can also be extended to rational functions in our vari¬ 
ables. An example for the outcome of our algorithm in the case of the non-selfadjoint 
rational function given by 


q(xi,x 2 ) := (0 
is shown in Figure [7] 


1 - \xx -ix 2 
1 ^ 


2> 1 



Figure 4. Brown measure (left) of p(x, y, z) = xyz — 2 yzx + zxy 
with x 7 y , z free semicircles, compared to histogram (right) of the 
complex eigenvalues of p(X, Y, Z) for independent Wigner matrices 
with N = 1000 



Figure 5. Brown measure (left) of p{x, y) = x + iy with x, y free 
semicircles, compared to histogram (right) of the complex eigen¬ 
values of p(X, Y) for independent Wigner matrices X and Y with 
N = 1000 
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Figure 6 . Brown measure (left) of p(x\, X2 , £3, X4 ) = X1X2 + 
X2X3 + X3X4 + X4X1 with xi, X2, X3, X4 free semicircles, compared to 
histogram (right) of the complex eigenvalues of p(Xi, X 2 , X 3 , X4) 
for independent Wigner matrices X-\, X - 2 , X 3 . X4 with N = 1000 



Figure 7. Brown measure of the non-selfadjoint polynomial 
q ( x \, X2 ) for free semicircular elements X \ and X2 
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